function traces = make_mogi_synthetic_waveform_fun(traces,...
    T, Q, A, step, noise, sourcenoise, onset,...
    mu, nu, lat, lon, elevation)
%create synthetic mogi waveforms
%(overrides data in existing iris trace structure, this is an easy way to
%make synthetic waveforms that have same instrument response and etc as
%desired real data)
% josh crozier 2020

g = 9.807; %gravity, leave positive

%get utc coords of source
zone = utmzone(lat,lon);
utmstruct = defaultm('utm');
utmstruct.zone = zone;
utmstruct.geoid = wgs84Ellipsoid;
utmstruct = defaultm(utmstruct);
[east, north, ~] = mfwdtran(utmstruct,lat,lon,elevation);

% separate by stations
stations = unique({traces.station});
keepind = true(size(traces));
for statin = 1:length(stations)
    coordin = cell(3,1);
    chanin = find(strcmp({traces.station},stations{statin}));
    for ii = 1:length(chanin)
        if strcmp(traces(chanin(ii)).channel,'HHE') %east
            coordin{1} = chanin(ii);
        elseif strcmp(traces(chanin(ii)).channel,'HHN') %north
            coordin{2} = chanin(ii);
        elseif strcmp(traces(chanin(ii)).channel,'HHZ') %vert
            coordin{3} = chanin(ii);
        else %delete trace
            keepind(chanin(ii)) = false;
        end
    end
%     coordin{1} = find(strcmp({traces.station},stations{statin}) & ...
%         strcmp({traces.channel},'HHE')); %east
%     coordin{2} = find(strcmp({traces.station},stations{statin}) & ...
%         strcmp({traces.channel},'HHN')); %north
%     coordin{3} = find(strcmp({traces.station},stations{statin}) & ...
%         strcmp({traces.channel},'HHZ')); %vertical
    
    %list of valid channel indices
    validin = false(3,1);
    if ~isempty(coordin{1})
        validin(1) = true;
    end
    if ~isempty(coordin{2})
        validin(2) = true;
    end
    if ~isempty(coordin{3})
        validin(3) = true;
    end
    
    %check that at least one valid channel present
    if any(validin)
        
        %make data even
        for cin = 1:length(coordin)
            if validin(cin)
                if mod(traces(coordin{cin}).sampleCount,2) == 1
                    traces(coordin{cin}).data = traces(coordin{cin}).data(1:end-1);
                    traces(coordin{cin}).sampleCount = length(traces(coordin{cin}).data);
                    traces(coordin{cin}).endTime = traces(coordin{cin}).startTime...
                        + traces(coordin{cin}).sampleCount/traces(coordin{cin}).sampleRate;
                end
            end
        end
    
        %select arbitrary channel
        tracein = coordin{find(validin,1)};
        
        %make synthetic source-time function
        time = seconds(1/traces(tracein).sampleRate*(0:traces(tracein).sampleCount-1));
        op = zeros(size(time));
        for Tin = 1:length(T)
            optemp = zeros(size(time));
            %sinusoid
            reltime = time(time > onset(Tin)) - onset(Tin);
            optemp(time > onset(Tin)) = ...
                A(Tin)*sin(reltime*2*pi/T(Tin));
            %decay
            reltime = time(time > onset(Tin) + T(Tin)/2) - onset(Tin) - T(Tin)/2;
            optemp(time > onset(Tin) + T(Tin)/2) = ...
                optemp(time > onset(Tin) + T(Tin)/2).*exp(-reltime*pi/(Q(Tin)*T(Tin)));
            
            %taper into smoother start (over half-wavelength of signal)
%             taperamp = 5*A(Tin)/(8*sqrt(2));
%             taperphase = -atan(4/3);
%             taperendin = find(time > onset(Tin) + T(Tin)/8, 1) - 1;
%             taperstartin = find(time > onset(Tin)...
%                 + (taperphase + pi/4)/(pi)*T(Tin), 1);
%             reltime = time(taperstartin:taperendin) - onset(Tin);
%             optemp(taperstartin:taperendin) = ...
%                 taperamp*sin(reltime*4*pi/T(Tin) + taperphase) + taperamp;
            
            %taper into smoother start (over full wavelength of signal)
            taperamp = A(Tin)/sqrt(2);
            taperphase = -pi/4;
            taperendin = find(time > onset(Tin) + T(Tin)/8, 1) - 1;
            taperstartin = find(time > onset(Tin)...
                + (taperphase)/(2*pi)*T(Tin), 1);
            reltime = time(taperstartin:taperendin) - onset(Tin);
            optemp(taperstartin:taperendin) = ...
                taperamp*sin(reltime*2*pi/T(Tin) + taperphase) + taperamp;
            
            %add step function (taper into start)
            taper = zeros(size(op));
            taperwave = T(Tin)/4;
            taperendin = find(time > onset(Tin) + taperwave/4, 1) - 1;
            taperstartin = find(time > onset(Tin) - taperwave/4, 1);
            reltime = time(taperstartin:taperendin) - onset(Tin);
            taper(taperstartin:taperendin) = ...
                step(Tin)/2*sin(reltime*2*pi/taperwave) + step(Tin)/2;
            taper(taperendin + 1:end) = step(Tin);
            optemp = optemp + taper;
            
            op = op+optemp;
        end
        
%         %add chirp
%         if ~isempty(chirpparam)
%             cp = chirp(tchirp,chirpparam{2},tchirp)
%         end
        
        %source noise
        op = op+max(A)*sourcenoise*randn(size(op));
        
        %highpass filter (prevents ULP noise)
%         op = highpass(op,1/(8*seconds(max(T))),traces(tracein).sampleRate,...
%             'Steepness',0.5);
        
        %get locations
        [stateast, statnorth, ~] = mfwdtran(utmstruct,...
            traces(tracein).latitude, traces(tracein).longitude,...
            traces(tracein).elevation);
        stateast_offset = stateast-east;
        statnorth_offset = statnorth-north;
        depth = traces(tracein).elevation - elevation;

        obs = zeros(3,1); %observation (station) locations
        obstilt = zeros(3,4);%observation locations for tilt
        obs(1,1) = stateast_offset;
        obs(2,1) = statnorth_offset;
        obs(3,1) = 0;
        obstilt(:,1) = obs(:,1) + [-1;0;0];
        obstilt(:,2) = obs(:,1) + [1;0;0];
        obstilt(:,3) = obs(:,1) + [0;-1;0];
        obstilt(:,4) = obs(:,1) + [0;1;0];

        %calculate greens functions for inflate-deflate motion
        mdl = [depth, 0, 0, 1]'; %[depth, east, north, volumechange]
        U = mogi_greensfunc(mdl,obs,mu,nu);
        Utilt = mogi_greensfunc(mdl,obstilt,mu,nu);
        G = [U(1,:).';U(2,:).';U(3,:).'];
        Gtilt = [atan((Utilt(3,2:4:end)-Utilt(3,1:4:end))/2).';...
            atan((Utilt(3,4:4:end)-Utilt(3,3:4:end))/2).';...
            0];

        %shift to fourier domain
        %interpolated one-sided freq (HZ)
        f = traces(tracein).sampleRate*(0:ceil(length(time)/2))/(length(time)).'; 
        Fop = fft(op);
        Fop = Fop(1:end/2+1);
        Fop(end) = conj(Fop(end));
        FU_pred = zeros(3,length(f));
%         FU_pred_tilt = zeros(size(FU_pred));

        %calculate displacements
        for fin = 2:length(f)
            M = [Fop(fin)];
            FU_pred(:,fin) = (G + g*Gtilt/(2*pi*1i*f(fin))^2)*M;
%             FU_pred_tilt(:,fin) = (g*Gtilt/(2*pi*1i*f(fin))^2)*M;
        end
%         FU_pred(:,1) = 0;
        
        %convert to velocity
%         'not converting to vel'
        for jj = 1:3
            FU_pred(jj,:) = FU_pred(jj,:).*(1i*f*2*pi);
        end

        %make 2-sided spectra
        FU_pred = [FU_pred(:,1:end-1),fliplr(conj(FU_pred(:,2:end)))];

        %Calculate U_pred
        U_pred = zeros(3,length(time));
        for pin = 1:3
            U_pred(pin,:) = real(ifft(FU_pred(pin,:)));
        end

        %add different noise to each station
        target_range = max([range(detrend(U_pred(1,:))),...
            range(detrend(U_pred(2,:))),...
            range(detrend(U_pred(3,:)))]);
        U_pred = U_pred + noise*target_range*randn(size(U_pred));

        %replace data with synthetics
        for cin = 1:length(coordin)
            if validin(cin)
                traces(coordin{cin}).data = U_pred(cin,:).';
                traces(coordin{cin}).sampleCount = length(traces(coordin{cin}).data).';
                traces(coordin{cin}).startTime = traces(tracein).startTime;
                traces(coordin{cin}).endTime = traces(tracein).endTime;
            end
        end
        
    end %end if any valid waveforms
end %end loop over waveforms

traces = traces(keepind); %delete bad traces

%convolve with instrument response
traces = add_response(traces);

end %end function
